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ABSTRACT:  Understanding  the  combustion  mechanism  of  jet  fuel  is  essential  for  the  design 
of  stable,  responsive,  and  highly  efficient  military  jet  engines.  However,  this  is  a  challenging 
mission  because  of  extreme  conditions  such  as  high  temperature,  supersonic  incoming  airflow, 
short  lifetime  of  intermediates,  the  complexity  of  combustion  reaction  networks,  and  the 
difficulty  to  identify  important  reaction  channels  in  combustion  chemistry.  In  the  present 
work,  a  physics-based  approach  (i.e.,  ab  initio  dynamics)  was  utilized  to  guide  the 
enumeration  and  characterization  of  chemical  reaction  networks  in  combustion,  minimize 
human  input,  and  avoid  on-the-fly  thermochemical  estimation.  Based  on  the  enumerated 
elementary  reactions,  a  framework  that  is  able  to  automatically  reduce  negligible  species  and 
elementary  reactions  was  suggested  and  hence  optimized  reaction  networks  were  generated. 
The  methodology  was  applied  to  study  methanol  oxidation,  prediction  of  ethylene  combustion 
mechanism  and  transition  states  of  different  types  of  chemical  reactions. 
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1.  Introduction 

Obtaining  more  powerful,  stable  and  responsive  jet  engines  is  always  an  essential 
need  for  the  nation.  One  important  component  of  meeting  this  need  is  for  researchers 
is  to  discover  and  understand  combustion  mechanism  of  jet  fuel  in  extreme  conditions 
such  as  supersonic  incoming  airflow  and  high  temperature.  This  is  a  challenging  task 
because  of  short  lifetime  of  intermediates,  the  complexity  of  combustion  reaction 
networks,  challenges  in  identifying  important  reaction  channels  in  combustion 
chemistry,  and  severely  lacking  experimental  methods  and  devices  that  are  able  to 
detect  and  monitor  a  significant  fraction  of  reactions  in  such  extreme  conditions. 

To  explore  combustion  mechanisms,  scientists  have  made  extensive  efforts  and 
great  progress.  For  example,  Smith  et  al.  provided  the  GRI  mechanism  for  natural  gas 
combustion  with  a  complete  set  of  elementary  reactions.1  The  researchers  of  the 
LLNL  combustion  projects  have  kept  producing  full  mechanisms  of  different 
hydrocarbon  fuels  for  over  thirty  years.2"4  Researchers  have  also  developed  different 
methods  and  algorithms  to  obtain  important  reaction  channels.  For  instance,  Broadbelt 
et  al.  suggested  a  matrix-based  reaction  network  generator  NETGEN  with  rate-based 
algorithm  to  optimize  reaction  channels.5"7  Green  and  co-workers  developed  Reaction 
Mechanism  Generator  (RMG)  to  produce  the  full  combustion  mechanisms  of  fuel, 
and  applied  rate-based  algorithm  and  sensitivity  analysis  to  predict  important  reaction 
pathways.8  Similarly,  there  are  combustion  reaction  network  generators  such  as 
COMGEN,  EXGAS,  and  KING  appearing  recently  as  well.9  All  these  reaction 
network  generators,  although  each  of  them  has  distinguishing  features,  share  common 
characteristics  and  generally  include  following  modules:  a  unique  topological 
representation  of  a  chemical  compound,  an  explicit  representation  of  reaction  rules,  a 
generation  scheme  through  the  application  of  reaction  rules,  and  a  gauge  to  optimize 
(reduce)  the  complex  reactions  networks.9  Human  inputs  like  reaction  rules  and  the 
method  for  optimization  of  reaction  networks  are  required  for  all  these  reaction 
network  generators.  In  addition,  to  predict  reaction  pathways,  on-the-fly  high-level 
quantum  calculations  are  often  needed  for  the  estimation  of  thermochemical  and  even 
kinetic  data.  Based  on  these  reasons  mentioned  above,  a  physics-based  approach  (i.e., 
ab  initio  dynamics)  was  suggested  to  guide  the  enumeration  and  characterization  of 
chemical  reaction  networks  in  combustion  and  minimize  human  inputs  and  avoid 
on-the-fly  thermochemical  estimation.  Furthermore,  a  framework  that  is  able  to 
automatically  reduce  negligible  species  and  elementary  reactions  was  suggested  and 
hence  optimized  reaction  networks  were  generated. 

In  the  present  work,  Car-Parrinello  molecular  dynamics  (CPMD)10  in  conjunction 
with  Metadynamics  (MetaD)  using  social  permutation  invariant  (SPRINT) 
coordinates  as  collective  variables,  as  the  physics-based  approach  mentioned  above, 
was  utilized  to  generate  combustion  network.  CPMD  is  a  planewave-based  ab  initio 
molecular  dynamics  method,  which  has  been  extensively  used  to  study  chemical 
reactions  because  of  its  faster  computational  capability  than  BOMD 
(Bom-Oppenheim  molecular  dynamics)11  and  reliability.12  To  overcome  the 
simulation  time  scale  problem  for  chemical  reactions,  the  enhanced  sampling  method 
MetaD  was  chosen  to  accelerate  the  exploration  on  potential  energy  surfaces  because 
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of  its  wide  and  successful  applications  to  chemical  reactions.1214  Meanwhile,  the 
SPRINT  coordinates  were  selected  as  generic  reaction  coordinates  because  it  was 
originated  from  graph  theory  and  includes  global  topological  information  of  all  atoms’ 
coordinates  in  the  system.  It  seems  to  have  ability  to  explore  all  reaction  channels  by 
solving  eigenvectors  of  contact  matrix  of  the  system.15 

By  using  the  biased  CPMD  approach  with  SPRINT  coordinates  this  AFOSR 
project  allowed  us  to  investigate  three  specific  research  areas:  methanol  oxidation, 
ethylene  combustion,  and  transition  state  searches  for  different  chemical  reactions. 
The  results  were  used  to  verify  the  practicability  of  the  suggested  physics-based 
approach.  And  hopefully  the  approach  can  be  extended  to  study  and  predict 
combustion  mechanism  of  jet  fuel  in  extreme  conditions  in  the  future. 

2.  Methodology  and  Computational  Details 

The  CPMD  v3.15.1  software  package10’16,  combined  with  PLUMED  vl.3,17  was 
used  to  run  all  Metadynamics18’19  molecular  simulations.  The  MetaD  technique  was 
selected  to  accelerate  sampling  of  rare  events.18  In  summary,  MetaD  forces  reactants 
to  react  within  the  underlying  phase  space  of  specific  degrees  of  freedom  by  adding  a 
history-dependent  Gaussian-shaped  bias  potential  on  coarse  descriptors  known  as 
collective  variables  (CVs).  Typically,  for  using  MetaD  to  study  chemical  reactions 
very  specific  CVs  related  to  the  researchers’  hypotheses  about  how  a  reaction  will 
happen  are  inputs  at  the  onset  of  the  simulation. 

The  choice  of  CVs  is  a  key  issue  in  MetaD  simulations.  Of  particular  importance 
for  generating  combustion  reaction  networks  is  the  requirement  that  the  CVs  must  be 
non-specific  because  what  kind  of  chemical  reactions  will  happen  are  not  predefined, 
and  can  be  driven  to  explore  all  elementary  chemical  reactions  and  reaction  channels 
by  adding  bias  on  them  to  explore  a  FES.  To  satisfy  this  requirement,  we  used  the 
SPRINT  coordinates,  which  include  the  information  of  all  atoms’  relative 
orientations  and  distances  to  define  the  topology  of  the  system  based  on  graph 
theory20'22,  as  the  CVs  in  CPMD  simulations. 

The  numbers  of  starting  configurations  of  interested  systems  were  set  large 
enough  so  that  sufficient  sampling  of  different  collision  pathways  were  able  to 
produce  major  reaction  channels.  All  elementary  reactions  and  channels  were 
summarized  and  generated  based  on  all  parallel  CPMD  simulations. 

3.  Discussion  and  Results 

To  monitor  the  time  evolution  of  various  species  during  combustion  and  verify  the 
capabilities  of  the  suggested  approach  to  enumerate  elementary  reactions,  a  series  of 
fuel  +  oxidation  molecules  simulations  were  set  up,  run,  and  investigated.  And  the 
results  were  summarized. 

First,  with  30  random  O2  configurations  relative  to  the  position  and  orientation  of 
CH3OH,  ICH3OH  +  IO2  simulations  were  performed  to  observe  the  methanol 
oxidation  process.  This  effectively  allowed  including  various  stochastic  effects  in  the 
system  since,  unlike  using  transition  state  theory,  the  outcomes  of  biased  CPMD 
simulations  potentially  could  be  dependent  on  starting  configurations  via  a  collision 
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trajectory. 

Figures  1  and  2  demonstrate  the  chemical  reactions  might  be  recorded  by 
monitoring  the  potential  energy  (the  Kohn-Sham  energy)  and  the  changes  in  the 
SPRINT  coordinates  of  involved  atoms.  Clearly,  the  first  energy  potential  energy 
barrier  is  for  the  bond  fission  of  C-0  in  methanol  as  shown  in  Figure  1,  and  ‘CPF 
radical  is  crucial  to  activate  the  stable  3C>2  molecule.  In  Figure  2,  around  time  9  ps,  all 
the  SPRINT  coordinates  are  in  a  narrow  range,  which  shows  that  the  reactants 
CH3OH  +  O2  are  close  to  each  other.  At  time  around  10.0  ps,  the  carbon  SPRINT 
coordinate  drops  first  due  to  the  C-0  bond  fission,  then  increase  by  1  because  of  the 
formation  of  CFpOO*  radical.  Similarly,  at  time  21.5  ps,  *H02  is  close  to  ‘CFFOH 
radical  and  CH2O  +  H2O2  are  formed,  therefore  their  SPRINT  coordinates  curves  are 
close  to  each  other  again.  In  fact,  we  can  observe  that  the  SPRINT  coordinates 
characteristically  change  along  with  chemical  reactions  in  all  30  simulations. 
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Figure  1.  The  time  evolution  of  elementary  reactions  and  Kohn-Sham  energy  of  one 

simulation  (#2)  at  temperature  1000K. 
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Figure  2.  Time  evolution  of  SPRINT  coordinates  of  one  simulation  (#2),  and  6  snapshots  at 
key  time  steps.  M-0  is  a  designation  for  the  two  atoms  initially  present  as  molecular  oxygen. 

To  investigate  the  properties  of  biased  CPMD  with  SPRINT  coordinates 
simulations  for  methanol  oxidation,  the  major  stable  species  were  tracked  during  the 
simulations.  As  shown  in  Figure  3,  the  counts  summarized  from  all  30  simulations  of 
selected  interesting  species  were  monitored.  The  motivation  for  the  aggregate 
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simulations  was  to  attempt  and  obtain  a  general  picture  about  what  the  MetaD 
trajectories  are  revealing  about  methanol  oxidation  given  that  the  biased  CVs  are 
totally  generic  about  the  types  of  transformations  that  could  happen.  Several  striking 
features  can  be  perceived  from  analysis  of  the  evolution  of  each  of  the  simulations. 
Oxidation  of  carbon-containing  species  is  gradually  observed,  and  in  a  production 
simulation  time  of  24  ps  all  fuel  molecules  completely  react  with  oxygen  in  all 
simulations.  As  shown  in  Figure  1,  following  methyl  radical,  the  methylperoxy 
(Ch^OO*)  are  the  largest  carbon-containing  intermediates  observed  and  play  crucial 
roles  for  succeeding  reactions.  In  the  later  stages  of  the  simulations,  a  mixture  of 
radicals  such  as  *CHO,  3,CH2,  ‘Ckf?,  *OH,  H02*  and  H*  appear.  Both  formaldehyde 
and  formyl  radical  are  produced  in  relatively  large  amounts,  suggesting  they  could  be 
a  gateway  between  the  substrate  and  final  products  CO2.  It  is  also  noteworthy  that 
hydroxyl  radical,  regarded  as  an  important  species  in  combustion  chemistry,  is  always 
existing  with  quite  stable  large  amount  (around  15)  in  the  production  runs.  In  the 
meantime,  the  numbers  of  H*  and  *CHO  steadily  increases  while  the  number  of  water 
regarded  as  a  quenching  molecule  in  combustion  keeps  growing.  Also  notable  is  that 
the  number  of  hydroperoxyl  (HCV)  radicals  is  quite  stable.15 
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Figure  3.  The  cumulative  time  evolution  of  count  numbers  of  selected  products  (H20  and  C02) 
and  radicals  with  largest  count  numbers  (*CHO,  3,CH2,  *CH3,  H*,  *OH,  H02*)  in  the  30 
CH30H+02  simulations  from  time  17.4  to  29.0  ps. 

To  further  test  the  practicability  of  the  new  method  in  combustion  science,  ethylene 
combustion  was  selected  as  the  target  since  there  are  still  some  open  questions  about 
its  mechanism  although  experimental  and  theoretical  scientists  have  extensively 
studied  it.  For  example,  Konnov  and  Xu  found  the  discrepancies  of  four  suggested 
ethylene  combustion  mechanisms.  By  studying  this  simplest  but  important  alkene 
with  the  new  approach,  the  re-discovered  mechanism  of  ethylene  combustion  was  to 
explain  the  conflicts  among  the  previous  four  mechanisms  for  same  fuel,  and 
therefore  verify  the  practicability  of  the  suggested  physics-based  approach. 

Three  model  systems  with  different  stoichiometry:  lethylene+102,  lethylene+202, 
and  2ethylene+102  were  set  up  for  CPMD  simulations.  The  temperature  was  set  at 
1000K.  By  combining  all  simulations  results  with  different  stoichiometry,  the  overall 
reaction  network  was  generated  as  shown  in  Figure  4.  The  total  connections  (reactions) 
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and  carbon-containing  species  are  386  and  116,  respectively.  Clearly,  even 
considering  the  most  fuel-rich  condition  (2ethylene+102),  it  was  confirmed  that  the 
C2  (2  carbon  atom  containing)  and  Cl  compounds  like  CH2CHO,  CH3CHO,  C2H5, 
CHO  and  CH3,  are  dominant  in  the  reaction  network  despite  of  so  many  C3  and  C4 
compounds  such  as  C3H6  and  C4H6  formed  in  ethylene  rich  conditions. 


Figure  4.  The  overall  graph  of  flux  of  carbon-containing  species  of  three  (1C2H4+1302, 
lC2H4+2302,  2C2H4+1302)  simulation  results  at  temperature  1000K.  The  area  and  color  of 
circles  represent  the  connectivity  of  different  carbon  species:  the  larger/darker  areas  have 
more  connections.  A  bidirectional  arrow  denotes  that  a  two-way  reaction  is  observed. 

Besides  the  combustion  chemistry  study,  because  of  the  unique  character  of  the  new 
approach,  i.e.  the  ability  to  identify  reaction  pathways  by  monitoring  the  potential 
energy  (EKs)  and  topological  descriptors  specific  reactants/products,  it  is  naturally  to 
apply  the  new  approach  to  detect  transition  states  for  different  chemical  reactions. 

A  well-studied  chemical  reaction  between  CH3CH2F  and  F'  in  gas  phase24,25  were 
investigated  using  the  new  approach.  There  are  6  possible  reaction  channels24  and  Sn2 
nucleophilic  substitution  CH3CH2F*  +  F'=  CH2CH2F  +  F*'  (*  denotes  the  substituted 
fluorine  atom)  reaction  is  one  of  them.  Without  specifying  which  atoms  were 
involved  in  the  reactions  (i.e.,  specifying  system-dependent  CVs),  20  biased  CPMD 
simulations  with  SPRINT  coordinates,  set  up  with  different  initial  configurations, 
were  run.  All  six  channels  were  successfully  reproduced  within  the  parallel 
simulations.  The  time  profiles  of  the  SPRINT  coordinates  and  Kohn-Sham  energy  of 
one  simulation  are  shown  in  Figure  5.  At  the  beginning  the  SPRINT  coordinate  of  the 
F  atom  (F2)  periodically  oscillates  between  3  and  4,  and  the  SPRINT  coordinate  of 
the  anion  F'  between  0  and  3;  at  time  about  1.05  ps  the  SPRINT  coordinates  of  the 
anion  fluorine  (aFl)  are  close  to  3  while  the  corresponding  Kohn-Sham  energy  has  a 
local  maximum  (the  black  circle).  The  transition  state  structure  for  Sn2  reaction  was 
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therefore  predicted  and  captured  as  shown  in  Figure  5.  The  captured  TS  structure 
from  CPMD  simulations  were  verified  by  directly  inputting  them  into  a  saddle  point 
search  in  Gaussian,  further  validated  by  confirming  the  presence  of  one  imaginary 
frequency,  and  coupled  with  IRC  calculations  to  verify  the  correct  reactants/products 
were  found.  Compared  to  Gaussian  2009  MP2/aug-cc-pVDZ  obtained  TS  structure, 
the  CPMD  obtained  geometry  agrees  well  with  MP2  results,  having  RMSD  0.02  A  for 
Sn2  TS  structures. 


Time(ps) 

Figure  5.  The  time  evolution  of  SPRINT  coordinates  and  Kohn-Sham  energy  of  a  Sn2 
reaction  CH3CH2F*  +  F'=  CH3CH2F  +  F*'at400K.  ‘aFl’  in  the  top  figure  denotes  the  anion 
F\  At  around  time  1.05  ps,  the  reaction  happens. 

Conclusions  and  Remarks 

Using  the  biased  ab  initio  method,  the  current  available  results  show  that  it  is 
possible  to  systematically  enumerate  all  elementary  reactions  in  the  simulations  of 
fuel  combustion,  generate  the  reaction  networks  with  carbon-containing  species,  and 
optimize  them  and  obtain  the  combustion  mechanisms.  For  example,  by  solving  the 
conflicts  among  the  previous  ethylene  combustion  mechanisms,  a  framework  that 
minimizes  human  input  and  maximizes  predictive  capabilities  has  been  created. 

By  correctly  identifying  transition  state  structure  of  CH3CH2F*  +  F'  =  CH2CH2F  + 
F*'  SN2  reaction,  it  shows  that  the  suggested  method  is  an  exciting  addition  to  the  list 
of  potential  uses  for  the  new  SPRINT  CVs,  and  paves  the  way  for  simultaneously 
exploring  both  the  interatomic  connectivity  in  complex  reacting  system  (i.e.,  a 
chemical  reaction  network)  along  with  accurate  representations  the  relevant  transition 
states  that  connect  stable  states  in  the  systems’  PESs.  By  applying  the  SPRINT 
coordinates  and  a  TS  identification  scheme  to  chemical  reactions,  the  new  approach 
seems  show  potential  to  recognize  TS  state  structure  accurately. 
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